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Preface 

This investigation determines an analytic solution to the 
force acting on a current loop over a non-uniform expanse of 
the geomagnetic field. Since the geomagnetic field contains 
energy and that energy is being used to make attitude 
adjustments to satellites, it should be possible to use that 
field energy as a sort of action-at-a-distance space 
thruster. Although the derivations for the force used a 
number of simplifying assumptions, the techniques developed 
are still applicable to a more detailed analytic solution for 
the force available through the interaction of a current loop 
and the geomagnetic field. 

I am deeply indebted to my thesis advisor, Lt Col Howard 
Evans, first, for the idea itself, and also, for his 
continued support and encouragement during the months of 
contemplation, analysis, evaluation, and writing that this 
effort represents. Without his superb instructional 
abilities and insights into the phenomenology of physics, 
this thesis would not have been completed. 

In addition, I would like to recognize my wife Olivia for 
her understanding and patience throughout my tenure at AFIT 
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Abstract 

This study determines analytic solutions to the 
interactive net force between the dipole created by a current 
loop and a non-uniform expanse of the geomagnetic field. A 
simple dipole approximation is used for the geomagnetic field 
and it is assumed that the current loop can be maintained in 
any desired orientation. Analytic solutions are found for 
four different orientations: the current loop parallel to 
the earth’s surface, the current loop perpendicular to the 
earth’s surface, the current loop dipole parallel to the 
instantaneous geomagnetic field, and an inertial orientation. 
A hypothetical satellite-current loop system is discussed and 
the action-at-a-distance magnetic interaction force is 
compared to aerodynamic drag and solar radiation pressure. 

The magnetic force is also compared to the force that some 
space thrusters produce. In terms of lifetime impulse or 
total impulse deliverable over one orbit, the magnetic force 
compares favorably to many thrusters in use today or under 
consideration. 
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GEOMAGNETIC FIELD ENERGY 
AS A SOURCE OF THRUST 


I. INTRODUCTION 


Background 

To accomplish their Department of Defense, scientific, or 
civil missions, satellites must periodically make attitude 
adjustments, compensate for orbital perturbations, or 
maneuver to different orbits. One constraint on their 
ability to accomplish these tasks, essentially a constraint 
on their lifetime, is the amount of fuel that can be carried 
on board. Systems currently in use include conventional 
chemical rockets, cold gas jets, and, in a few situations, 
electric propulsion. Once the fuel supply for these 
activities is exhausted, the satellites cease to be useful 
and become expensive orbiting junk. An energy source is 
needed that does not require the carriage of an on-board fuel 
supply. 

The earth is surrounded by a magnetic field; this field 
contains energy (B 2 /2yo). By taking advantage of this 
energy, the amount of on-board fuel consumed by a satellite 
can be reduced and the lifetime increased. While magnetic 
torques have been used to adjust satellite attitude (Eller, 
1983:315; Kaplan, 1976:196), there has been no attempt to use 
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the energy contained in the geomagnetic field to change the 
satellite's orbital elements. Regardless of the means used 
in moving a satellite, all require an exchange of momentum 
between an energy source and the space vehicle. 
Traditionally, any changes to a satellite’s orbit have 
required the use of thrusters and a limited supply of 
reaction mass. Since no reaction mass will be used in this 
proposed system, the momentum exchange will take place 
between the earth and the satellite through the interaction 
of their respective magnetic fields. 

Statement of Purpose 

This thesis will analyze the interaction between a 
magnetic dipole field generated by a current loop around an 
artificial satellite and a non-uniform expanse of the 
geomagnetic field. Specifically, it is to quantify the net 
force available through the interaction of the two fields. 
This action-at-a-distance force will then be compared to two 
orbit degradation perturbative forces: atmospheric drag and 
solar radiation pressure. A comparison will also be made to 
representative space thrusters currently in use or under 
consideration. In addition, the force will be evaluated in 
terms of its potential capability to make minor orbit 
alterations. 
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Scppq 


It is known that satellite attitude and spin rate can be 
changed by the earth’s magnetic field as a result of the 
torque a magnetic dipole experiences in a uniform magnetic 
field. This thesis will look at the net force available to 
act on a satellite in a non-uniform external magnetic field 
and the possibility of maintaining an artificial satellite’s 
orbital characteristics using geomagnetic field energy in the 
face of natural decay process or to make minor adjustments to 
the orbital elements of an artificial satellite. 

This research will entail a determination or analysis of: 

-The characteristics of the earth’s magnetic field 
-The interaction between a magnetic dipole and a 
uniform magnetic field 

-The interaction between a magnetic dipole and a non- 
uniform magnetic field 

—Specifically the net force on a satellite - as 
a result of a large rigid current loop about 
the satellite - in the earth’s magnetic field 
-The effect on a satellite’s orbit by the application 
of such a force 

To achieve this, certain assumptions and limitations must 
be established: 

-The system starts out in a stable orbit 

—Only polar orbits will be considered in the 
derivations 
-A non-rotating earth 

-The geomagnetic field changes at a predictable rate 
-The current carrying loop about the satellite is 
rigid 

-Power requirements can be met 

--No specific hardware will be designed 
-The current carrying loop can always be maintained 
in a given orientation relative to the geomagnetic 
field 
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Scope 


A literature search and review are necessary to resolve 
several preliminary and fundamental knowledge requirements 
for this thesis. The four primary categories of search and 
review are orbital mechanics, magnetostatics, the structure 
of the geomagnetic field, and present usages of the 
geomagnetic field. 

The topics of orbital mechanics and magnetostatics are 
well developed and thoroughly covered in textbooks. The 
specific area of interest in orbital mechanics is the change 
that takes place in a satellite’s orbit by the continuous 
application of a very small force—rather than the single 
impulse approximation used with conventional chemical 
rockets. Small, non impulsive forces are generally dealt 
with through perturbative techniques. One such technique, 
with application to this thesis and presented in Appendix A, 
is Euler’s Variation of Parameters. 

Within magnetostatics, several preliminary problems need 
to be understood prior to attempting the computations for 
this thesis. These are: 

-The magnetic dipole generated by a current loop 

-The interaction between a magnetic dipole and a 
unifbCB magnetic field 

-The interaction between a magnetic dipole and a non- 
uniform magnetic field 
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The first is covered in most undergraduate textbooks on 
general physics. The second is covered in textbooks dealing 
with electromagnetic theory and represents the mathematical 
concept behind the use of magnetic torques for satellite 
attitude adjustment (Kaplan, 1976:196f). The last problem is 
not often addressed but when addressed it is in terms of 
microscopic particles. 

The third category of search and review is the requirement 
to find an acceptable model of the geomagnetic field. This 
model will be integrated into the problem of the interaction 
between a magnetic dipole and a non-uniform magnetic field. 

The final area of review is the current use of the 
geomagnetic field for satellite attitude adjustment. Several 
related areas that will not be investigated due to 
simplifying assumptions, but would be necessary for a 
practical engineering design resolution to the thesis 
question, include satellite power supply, structural design 
requirements, satellite attitude stabilization schemes, and 
launcher constraints. 

Orbital Mechanics 

Since the action-at-a-distance force that may be derived 
from the interaction of the satellite magnetic dipole field 
and the geomagnetic field will be small, the application of 
the force must be maintained for finite periods and 
preferably in the most intense region of the geomagnetic 
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field as possible. This non-impulsive application of force 
suggests a similarity to low-thrust systems such as ion 
propulsion (Stuhlinger, 1964:84,89). Low thrust problems 
must be analyzed using a perturbation approach such as 
Cowell's, Encke’s, or Variation of Parameters (Bate and 
others, 1971:425). While the optimum impulsive thrust 
transfer between two coplanar circular orbits is a Hohmann 
transfer, the optimum infinitesimal thrust transfer is a 
spiral. Close to the earth the gravitational acceleration is 
much larger than the thrust acceleration and the spiral will 
be very tight. A graphic comparison of a Hohmann transfer 
and a low-thrust transfer is shown in Figure 1. 



Figure 1. Comparison of a Hohmann Transfer Between Circular 
Orbits (a) and a Low-Thrust Spiral Transfer Using 
Infinitesimal Thrust (b) 

(Adapted from Edelbaum:117) 
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Over time the acceleration effects accumulate and the 
orbit expands more quickly. A general assumption in low- 
thrust transfer maneuvers is an acceleration that is 
approximately constant in magnitude and tangential in 
direction. A review of known solutions to some trajectory 
optimization problems using impulsive thrust shows that the 
same orbital transfers are possible using infinitesimal 
thrust (Edelbaum, 1966:113). 

M.ajsmg.fcg.S tat i c§ 

Unlike electrostatics, which can point to a unit charge as 
its fundamental building block, magnetostatics has no 
magnetic "monopole" that fulfills the same role (Cheston, 
1964:149). The fundamental magnetic material is the "motion 
of charges" (Dart, 1966:66). The motion of electric charges 
produces a magnetic field. The magnetic field intensity-- 
often referred to as magnetic induction—is a vector field 
having both direction and magnitude. It is proportional to 
the current and inversely proportional to the distance from 
the current (Dart, 1966:68). 

When a current is introduced into a magnetic field, there 
is an interaction. This interaction is felt as a force at 
right angles to both the current and the field. The 
magnitude of the force is the product of the external 
magnetic field intensity and the current (Halliday A Resnick, 
1970:559). Figure 2 shows the orientation between the 
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magnetic field, current, and direction of the force, using 
the right-hand rule. 

If the current forms a circular loop, the magnetic field 
generated is the same field that would be generated by a 
magnetic dipole oriented perpendicular to the plane of the 
circular loop (Halliday & Resnick, 1970:571). If the 
current-carrying loop is introduced into an external uniform 
magnetic field, it will experience zero net force and hence 
no displacement. However, if the dipole moment of the loop 
is at an angle (8) to the direction of the external magnetic 
field, there is a net torque; and the loop will rotate if 
free to do so as shown in Figure 3. (Halliday & Resnick, 
1970:543). 



Figure 2. Relative Direction of Magnetic Field, Current & 
Force 









Figure 3. Forces on a Current Loop in a Uniform Magnetic 
Field (Adapted from Halliday & Resnick:542) 


The limited discussion encountered in the literature of a 
magnetic dipole in a non-uniform magnetic field provided a 
generalized equation for the force on a localized current 
distribution in a slowly varying external magnetic field 
(Jackson, 1967:167). The author’s discussion applied the 
equation to a charged particle—such as might be found in the 
Van Allen belts surrounding the earth--and provided the basis 
for the magnetic mirrors discussed in Evans (1985:109). A 
schematic of the forces acting on a magnetic dipole in a 
uniform and non-uniform magnetic field are shown in Figure 4. 
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Figure 4. Forces on the Effective Dipole of a Current Loop in 
a Uniform (a) and a Non-Uniform (b) Magnetic Field 


Geomagnetic Models 

A magnetic field in free space is usually given as a 
vector potential; however, it can also be represented in 
other formats: scalar potentials, orthogonal vectors, Euler 
potentials or local expansions of the magnetic field around a 
reference point. All of these additional forms, however, can 
be converted back to the vector potential format (Stern, 

19 7 5 : 2 f). 

The geomagnetic field, like any magnetic field, is a 
vector field, and magnetic field computations can be based on 
the additive properties of vectors (Klump&r, 1982:361). 
Ideally, the magnetic field at any given point on the surface 
of the earth or at a satellite’s orbit could be calculated by 
summing up all the contributions from the currents that 









influence that point. However, since the earth and its near 
environment are not a homogeneous mass nor a simple linear 
current distribution system, the actual vector field is 
complex and can only be approximated. A common method of 
approximating a complex function is to use an infinite series 
expansion. The series can be expanded as necessary to find 
the dominate terms or to reach a required level of precision. 
Representations of the geomagnetic field potential are 
presented, in most cases, as truncated infinite series. The 
International Geomagnetic Reference Field 1975 (IGRF 1975.0) 
is used to facilitate comparisons of different models. 
Represented by a truncated harmonic series of 80 terms, it 
was constructed by averaging the coefficients of several 
models. 

The simplest representation of the geomagnetic field is 
the field created by the first term of the expansion. 

Referred to as the axial dipole model, it is the field 
created by a dipole whose axis is the earth’s rotation axis. 
This model can be further simplified by the assumption of a 
perfectly circular orbit (i.e. a perfectly spherical earth). 

If the expansion is carried out to three terras, the model 
now represents the potential of a dipole whose axis is tilted 
11.4 degrees to the earth’s axis of rotation. This is the 
oblique or tilted dipole. Further enhancement of the oblique 
dipole can account for the center of the earth’s dipole and 
the earth’s center of mass being separated by up to 500 
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kilometers (Evans, 1985:93). Various sources and models give 
displacement distances of 300, 342, 467, and 500 kilometers 
(Stacy, 1969:129; Alfven, 1963:5; Wertz, 1985:779; Evans, 
1985:93). As the radial distance from the earth increases, 
the dipole approximation becomes increasingly representative 
of the actual field. However, beyond about two and a half 
earth radii the dipole nature of the geomagnetic field is 
distorted by solar activity and other phenomena (Eller, 
1983:317) . 

Models of magnetic fields based on the mathematics of 
spherical harmonics [related to orthogonal vectors] take the 
magnetic field at the surface of the earth and expand it 
radially (Stern, 1975:6). Again, beyond two or three earth 
radii, the solar wind begins to distort the geomagnetic field 
and this technique has decreasing utility (Eller, 1983:316f). 

Quantitative models of the geomagnetic field tend to fall 
into two classes: internal models, where the field is due to 
sources inside the earth and is derived from some type of 
direct measurements (spherical harmonics and the axial dipole 
model); and external models that attempt to explain the 
field by adding the effects of currents above the earth’s 
surface (Mead A Fairfield, 1975:523). As data has become 
available, models have been--and probably will continue to 
be--made based on actual satellite measurements of these 
external currents. 
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One external model often referenced, and the first to be 
directly based on field measurements, is the Mead-Fairfield 
1972 model. This model was derived from a least-squares fit 
to 12,616 vector field averages. The measurements were made 
by four Explorer satellites between 1966 and 1972. The model 
provides four sets of coefficients to account for 
"superquiet, quiet, disturbed, and superdisturbed" (Mead & 
Fairfield, 1975:526) geomagnetic activity. It is considered 
applicable from five to seventeen earth radii (Fairfield & 
Mead, 1975:535). 

In 1979 a global vector survey of the geomagnetic field 
was initiated with the launch of MAGSAT. Its mission was to 
measure the main geomagnetic field as well as the variations 
caused by crustal anomalies. Due to the sensitivity of the 
on-board instrumentation, it was also able to investigate 
other external currents such as the equatorial, auroral zone 
and polar cap, and high latitude ionosphere and magnetosphere 
currents. From these measurements a model was developed to 
account for geomagnetic field perturbations due to various 
currents in the ionosphere and magnetosphere (Klumpar, 
1982:363). 

Actual Usages o f the Geomagnetic Field 

The resultant torque mentioned previously of a magnetic 
dipole in a uniform magnetic field is the principle by which 
magnetic torques are employed to adjust satellite attitude. 
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In the mid-1960’s Harold Perkel conceived of a system to 
control transverse momentum components using the interaction 
of the geomagnetic field and magnetic torquers on board a 
satellite. When a magnetic coil’s dipole is aligned parallel 
to the spin axis of a satellite, precession control torques 
are generated through interaction with the geomagnetic field. 
By switching the sign of the magnetic coil’s magnetic moment 
at appropriate points in the orbit, the satellite’s spin axis 
can be precessed in any desired direction (Herman, 1977:5). 
The concept was originally used in the TIROS weather 
satellite program and later incorporated into the synchronous 
RCA communications satellites (Kaplan, 1976:196). 

This tendency for a magnetic dipole to align itself with 
an external magnetic field is also one of the primary sources 
of disturbance torques on satellites (Eller, 1983:316). 
Vanguard 1 despun from 2.7 to 0.2 revolutions per second in 
two years due to interaction between on-board magnetic fields 
and the geomagnetic field (NASA, 1969:2). The Navigation 
Demonstration Satellite (NDS-2) built by Rockwell for the 
Global Positioning System (GPS) was originally built to dump 
momentum with hydrazine thrusters that torque the satellite 
in the proper direction. However, "Since 1980 the satellites 
of the GPS system have dumped momentum magnetically, thus 
avoiding ephemeris disturbance and hydrazine consumption" 
(Eller, 1983:316). This represents an example of real fuel 
savings through the use of the geomagnetic field energy. 
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The calculation of net force acting on a magnetic dipole 


in a non-uniform magnetic field is representative of the 
problem of a satellite and its magnetic field operating in 
earth orbit. The simplest representations of the geomagnetic 
field are the axial dipole models. These are considered 
accurate to about three earth radii. For proof of concept 
and order of magnitude calculations in low to medium earth 
orbit, a first order spherical harmonic model is a valid 
approach. For the purposes of this thesis, therefore, an 
axial model will be used since only low to medium orbits will 
be considered and the geomagnetic field intensity diminishes 
as the cube of the distance away from the earth. 
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This section will introduce and develop the physics and 
mathematical equations behind magnetic dipoles and fields and 
their interaction. Following this discussion will be a 
restatement of the two-body problem and a brief look at 
orbital mechanics. The closing section will discuss the 
equations for geomagnetic field models that exist and are 
pertinent to this thesis. 

Magnetostatics 

A starting point for this study is an investigation into 
the characteristics of magnetic fields. A magnetic field 
results from the motion of electrically charged particles 
(such as an electric current in a wire as shown in Figure 5) 
or the presence of a magnet, and at a given point generates a 
sideways deflecting force on another electrically charged 
particle that has an arbitrary initial direction and speed 
(vector) through that point. Figure 6 shows the curved paths 
positive and negative electrons take in the presence of a 
magnetic field as a result of the sideways deflecting force. 

The most familiar manifestation of a magnetic field is the 
torque (twisting or rotating motion) that a small magnet, 
such as a compass needle, experiences when placed in such a 
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region. The compass needle will align itself in a north- 
south direction in the earth’s magnetic field or at right 
angles to a current carrying wire. 



Figure 5. Magnetic Field Around a Current Carrying Wire 
(Reprinted from Halliday and Resnick:561) 



Figure 6. Paths of Positive and Negative Electrons Deflected 
in Opposite Direction in a Magnetic Field 
(Reprinted from Halliday and Resnick:540) 
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The torque is proportional to the field intensity and to 
the sine of the angle between the axis of the magnet (the 
line joining the poles) and the field direction. The 
direction is determined by the equilibrium position that the 
north pole of a compass needle will reach when placed in the 
magnetic field. 

Unlike the unit of electrical charge that exists, there is 
no known magnetic monopole. The simplest magnetic structure 
is the magnetic dipole, characterized by magnetic dipole 
moment (p). The simplest representation is the familiar bar 
magnet, but magnetic dipoles are also generated by current 
loops and solenoids. 

A mathematical relationship between the magnetic field 
vector B and the current i is given by Amperes Law. 

/fi • dL = po i ( 1 ) 

where 

B = magnetic field evaluated around a closed path in 
space (tesla) 

dL = displacement vector around a closed path (meters) 

i = current contained within the path of integration 
(amperes) 

po = permeability constant 

Its usefulness is limited in that the symmetry of the current 
distribution must be high enough to permit evaluation of the 
integral such as represented by Figure 7a. Even moderately 
complex paths can be difficult to integrate (Figure 7b). 
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Figure 7. Paths of Integration in a Magnetic Field 
(Reprinted from Halliday and Resnick:560) 


An example of its application can be found in determining 
the magnetic field near a long wire. 

B = po i/2itr (2 ) 

where r is the radial distance from the wire. 

To compute fi at any point due to an arbitrary current 
distribution the Biot-Savart Law states 

dfi = (po i/4n ) [ {dh x r)/r 3 )] (3a) 

or 

dB * (poi/4re)[dLsin(0)/r 2 ] (3b) 

so that 

fi =/dfi (4) 

In Figure 8, "P" is a point where dfi is the magnetic field 
associated with current element dfi. The magnitude of dfi is 
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given by Eq (3b) and £ and 6 are the displacement vector and 
direction to "P" from dL* 



Figure 8. Schematic representation of the Biot-Savart Law 
where Current Element d]* Establishes a Magnetic 
Field Contribution dfi at "P" 

(Adapted from Halliday and Resnick:569) 


Therefore a current distribution such as a current loop, 
generates a magnetic field—the same field that a magnetic 
dipole aligned perpendicular to the plane of the loop would 
generate. Application of the Biot-Savart Law gives the 
magnetic field along the axis of a circular current loop to 
be 

B = poiR 2 /2(R 2 + x 2 ) 3 /* ( 5 ) 

where 

R = radius of the loop 

x = distance from the center of the loop to a point along 
the axis as shown in Figure 9 
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if x >> R 


B = poiR 2 /2x 3 


( 6 ) 



Figure 9. Magnetic Field Along the Axis of a Current Loop 
(Reprinted from Halliday and Resnick:570) 

Remembering that itR 2 is the area ’A’ of the current loop 
and defining the magnetic dipole (p) as the current in the 
loop (i) times the area of the loop (A) times the number of 
turns of wire in the loop (N) yields: 

B = poNiA/2nx 3 = pop/2xx 3 (7) 

Since the geomagnetic field may be approximated by that of 
a simple magnetic dipole field, the magnetic field around a 
magnetic dipole is of interest. The equations describing a 
magnetic dipole field are derived from the magnetic potential 
in the vicinity of a magnetic dipole. 

Referring to Figure 10, let P be a point within the 
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influence of a simple dipole located at the origin, with P’s 
location described in spherical coordinates (r,0,0). 



Figure 10. Geometry of Potential Around a Magnetic Dipole 

The magnetic moment (y) is given by the separation ly) 
between the poles times their pole strength (± m). 

y = my (8) 

The potential U at P is then: 

U = (-m/ri) + (m/r 2 ) (9a) 

= -(m)/([r 2 +(y/2) 2 -rycost0 ) ] 1 / 2 } 

+ (m)/{[r*+(y/2) 2 +rycos(0)] 1 / 2 } (9b) 

If the separation between the poles is much less then the 
distance to P, that is y<<r, the potential can be given by: 

U = -mycos(0)/r 2 = -ycos(G)/r 2 (10) 

The magnetic field at P is then: 







Letting X = (n/2)-6 be the magnetic latitude, the magnetic 
field can be defined in terms of its local horizontal (B ) 
and vertical (Br) components as: 

B X = (l/r)OU/38) = (2p/r 3 )(cos(\)/2) (12) 

Br= 3U/Br = 2psin(X)/r 3 (13) 

The total field strength is then: 

B = (B\ 2 +Br 2 ) 1 / 2 = (p/r 3 )[l + 3sin 2 (X)] 1 / 2 (14) 

Now consider the affect magnetic fields have on objects 
within their influence. A positive test charge qo moving 
with velocity v in the presence of a magnetic field B 
experiences a force given by: 

F = qo v x B (15) 

where 

F = force (newtons) 

v a velocity of charged particle (meters/second) 
qo = charge (coulombs) 

B = magnetic field (tesla) 

The units for B stem from this relationship: 

(nt/coul)/(meter/sec) = weber/meter 2 = tesla 

1 tesla = 1 weber/meter 2 = 10000 gauss 

If a charged particle is moving through a region in which 

both a magnetic field and an electric field exist, the force 

experienced is given by the Lorentz relation: 

E = qof + qov x B (16) 


The fact that the magnetic force is always at right angles 
to the direction of motion means that (for steady magnetic 
fields) the work done by this force on the particle is 
zero. ... Thus a steady magnetic field cannot change the 
kinetic energy of a moving charge; it can only deflect it 
sideways. (Halliday & Resnick, 1970:539). 
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Therefore it is reasonable to assume that since a magnetic 
field exerts a sideways deflecting force on a charged 
particle it would likewise exert a sideways force on a 
current carrying wire. This force is given by: 

E = Hi x fc (17) 


where 

I* = displacement vector pointing in the current 
direction (meters) 
i = current (amperes) 

£ * magnetic field intensity (tesla) 

If a current carrying loop is inserted in a uniform 

magnetic field it will experience a torque (t) given by: 

t = a x 1 ( 18 ) 


with magnitude 


t = pBsin(9) 


( 19) 


with 

U = magnitude of magnetic moment of loop = iA 
9 = angle between vectors 14 and £ 

The torque will act on every turn of a multiply coiled loop 

so that 


t a NiABsin(Q) 


( 20 ) 


where 

N = number of coils in the loop 
i = current (amperes) 

A = area of the loop 

9 = angle between the magnetic field and the magnetic 
dipole moment of the loop 
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Two Body Problem 


In the two body problem of an earth satellite orbiting 
about the much larger mass of the earth, the gravitational 
attraction (Fg) between the two bodies is given by: 

Fg = -GMm/r 2 (21) 

where 

M = mass of the earth (kg) 

m = mass of the satellite (kg) 

r = distance of satellite from center of the earth (m) 

G = gravitational constant (Nt*m/kg 2 ) 

This is based on the following two assumptions: 

1) The earth and satellite are spherically symmetric 
so they can be considered as point masses. 

2) There are no external nor internal forces acting 
on the pair aside from the gravitational force 
between them. 

If these assumptions are allowed, as well as assuming the 
satellite is already in a stable orbit, its orbit is governed 
by the following equations of motion: 

X = -G(M+m)x/r 3 (22) 

Since, in our system M>>m , it is convenient to say G(M+m) = 
GM and define GM = k, then 

X + kx/r 3 =0 (23) 

A satellite’s position can, in the appropriate coordinate 
system be defined by a vector x and its velocity by the 
vector v. The moment of momentum (angular momentum per unit 
mass) vector (h) is defined by h = £ x v and for a given 
orbit is constant. 
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Following the derivation in Thomson and defining u = 1/r 
the differential equation for the orbit becomes 

d 2 u/d0 2 + u = k/h 2 (24) 

This is a second order differential equation requiring two 
arbitrary constants in the general solution. 

u = k/h 2 + C cos (0 -0o) (25) 

With the kinetic energy per unit mass given by v 2 /2 and the 
potential energy per unit mass by -k/r the total energy per 
unit mass (E) is: 

E = (v 2 /2 ) - (k/r) (26) 

or 

E = -k/2a (27) 

where 

a = the semi-major axis of a conic section 
Choosing 0o = 0 in Eq (25), the solution to the orbit equation 
can be written 

u = k/h 2 (1 + e cos(0)] (28) 

where e is the eccentricity for any conic orbit and is given 
by: 

e = (1 + 2Eh 2 /k 2 ) l /2 (29) 

Using the above defined characteristics any orbit can be 
fully described in terms of shape, size and orientation by 
five independent quantities referred to as "orbital 
elements". 
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These are: 


1) e = eccentricity 

2) a = semi-major axis 

3) i = inclination 

4) x = longitude of ascending node 

5) w = argument of periapsis 

If the satellite position must be pinpointed at a particular 
time a sixth element is also needed. 

6) T = time of periapsis passage 

Other elements are sometimes used such as the semi-latus 
rectum (p) which is derived from the semi-major axis and the 
eccentricity, or the longitude of periapsis which is the 
longitude of the ascending node plus the argument of 
periapsis. These relationships are shown graphically in 
Figure 11. 

The preceding discussion is dependent upon the two 
assumptions previously mentioned. To account for the 
perturbing effects of other masses (sun, moon, planets, etc.) 
a series of terms must be added to Eq (23) so that we have: 

X+( k£/r 3 )+E n j = 3 Gmj [ (xi 2 /rj 2 3 ) - (rJ i/rj l 3 ) 1 =0 (30) 

The external accelerations caused by these other masses 
are small relative to the acceleration by the earth’s 
gravitational field (Bate and others, 1971:11) and for this 
thesis will be ignored. Other interactions, not considered 
Keplerian, include perturbations due to the nonspherical 
earth, atmospheric drag, solar radiation forces, solar wind, 
eddy currents in the satellite, electromagnetic induction, 
cosmic dust etc. The interaction of a controlled, internal 
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current with the earth’s magnetic field is the specific 
perturbation to a satellite’s orbit to be considered in this 
thesis. The variation in the elements that define the orbit 
over time will be determined based solely on the magnetic 
field interaction. 

Techniques to account for these perturbations include 
Variations of Parameters or Elements, Cowells Method, and 
Encke’s Method . This thesis will use the Variation of 
Parameters or Elements as discussed in Bate (1971:396) to 
outline the procedures to determine the effects of the 
magnetic field interaction on a satellite’s orbital elements. 
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Figure 11 


Orbital Elements 

(Reprinted from Bate and others:59) 
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The exact cause of the earth’s magnetic field is not 
known. It is believed to be a result of currents within the 
earth’s core or mantle; the effect of currents located above 
the surface such as in the ionosphere, or a combination 
thereof. Gilbert proposed that the earth was a great 
spherical lodestone. This was discounted due to the lack of 
lodestone located near the surface. Later Gauss introduced 
the spherical harmonic expansion of a scalar potential 
analysis that fit the measurable portions of the geomagnetic 
field. 

Since the field rotates with the earth it is reasonable to 
assume that most of the field originates from some internal 
phenomena such as magnetization or a dynamo effect. A field 
internal to the earth can be described by a solution to a 
boundary value problem. There are no surface currents on the 
earth so curl £ = 0 [note: Klumpar (1982:361) points out 
there are external sources, and Mead & Fairfield conclude 
that curl £ ^ 0 --though by a very small quantity, but 
Stern (1975:7) points out that internal sources account for 
99% of the field]. This lets the field be expressed as the 
gradient of a scalar potential g = -grad U . Basic 
electromagnetic theory says that with no known magnetic 
monopole, we also have div £ = 0 . Combining the last two 

equations yields Laplace’s equation: 

curl grad U = 0 (31) 
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The solutions of this equation are well known and can be 
expressed in terms of functions known as spherical harmonics 
which are given in units of geocentric coordinates as 

U(R,e,t)=aI (a/r)"* 1 [ g*n cos( m# )+h*n sin(m<t) ]P"n ( Q ) (32) 

n=l m=0 

where 

a = equatorial radius 

g“n, h"n = Gaussian coefficients (determined 

empirically by least squares fit of 
measured values) 

(R,d,$) = geocentric distance, coelevation, east longitude 
from Greenwich which defines a point in space 
P* n( 9 ) = associated Legendre functions 

When n=l the expansion yields a simple dipole, with n=2 
the expansion represents a quadrupole, etc. Surface 
indications such as the Brazilian anomaly (Wertz, 1985:115) 
indicate the dipole is offset from the earth’s center by 
approximately 474 kilometers (Chapman and Bartels, 1940:12). 
Therefore the eccentric nature of the dipole is accounted for 
by the quadrupole expansion. With n=3 the field is that 
given by an octopole. The IGRF model is a result of n=m=8. 

The components of the field are given by the respective 
partial derivatives. 

Various analytical attempts have been made to calculate 
the coefficients but with little success. Field measurements 
have yielded coefficients empirically through least squares 
fits. While accurate for the data from which they were 
calculated, the earth’s magnetic pole strength and 
orientation are subject to secular drift. Consequently the 
accuracy of the models degrades over time. 
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Since the dipole field decreases as 1/R 3 and the 
quadrupole decreases by 1/R 4 and higher multipoles decrease 
even more rapidly, expansions to fewer terms are justified 
for some calculations. Wertz (1985:119) shows the effects of 
field truncation and angular errors. As altitude increases 
the appearance of the geomagnetic field more and more closely 
resembles that of a magnetic dipole. This thesis is a 
conceptual analysis attempting to quantify the net force 
available through the interaction of a satellite magnetic 
dipole and the geomagnetic field. Instantaneous precise 
forces are not germane and would only be as accurate as the 
currency of the coefficients used and will therefore not be 
attempted. The geomagnetic field will be approximated by an 
axial dipole model. 
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IV. DEVELOPMENT 


The force that the geomagnetic field exerts on a large 
orbiting current loop will vary depending on the current the 
loop is carrying, its radius, and position and orientation 
relative to the geomagnetic field. A central assumption of 
this thesis is that the current loop can be maintained in any 
desired orientation relative to the earth and the geomagnetic 
field throughout its orbit. 

Four orientations will be discussed and analyzed. Case 1 
will align the satellite’s magnetic dipole parallel to the 
radial component of the geomagnetic field, so tha*- the plane 
of the current loop is parallel to the earth’s surface. Case 
II will look at the force as a result of an alignment of the 
current loop’s dipole with the theta component of the 
geomagnetic field, so that the current loop is always 
perpendicular to the earth’s surface. Case III will 
investigate the force as a result of alignment of the current 
loop’s dipole with the instantaneous geomagnetic field at 
each point in the satellite’s orbit. In Case IV the 
satellite will maintain a constant inertial orientation. 

This fourth orientation can be generalized to any arbitrary 
orientation. Figure 12 shows the four cases and their 
orientation with respect to the earth, each other and the 
components of the geomagnetic field. 
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The equations representing the force on the current loop 
will now be derived as functions of the radial distance from 
the center of the earth and the co-latitude position. While 
the geomagnetic field is strongest and a greater gradient 
exists near the poles, the net effective force is not 
necessarily greatest at the poles. However, since the force 
will act in the theta or radial direction the interactive 
force is most advantageously applied to a satellite in a 
polar orbit. 

The force on each section of the current loop is due to 
the interaction between the current displacement vector for 
that section and the geomagnetic field vector acting on that 
section. The direction of the force is perpendicular to both 
the current vector and geomagnetic field. This implies that 
the radially directed component of the force ( Fr ) is a result 
of the theta component of the geomagnetic field (Bt). 
Likewise, the force component in the theta direction ( Ft ) is 
due to the radial component of the geomagnetic field ( Br). 

The force in the phi (longitude) direction (Fp ) will be shown 
to be zero. 

It is the force on the loop itself that is of interest, 
but the position of the loop will be defined by the location 
of its center only, that is, a radial distance (Ro) and co¬ 
latitude theta (9o). The geomagnetic field does not change 
with the longitude, therefore the longitudinal position of 
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the center (4o) does not need to be considered in this 
derivation since it does not contribute to the solution. The 
position of the center of the loop (Ro,0o,4o) will be referred 
to by position vector fio. The magnetic field interacting 
with the loop at location £ is the magnetic field at the 
radius, r, from the center of the loop. This is given by the 
gradient of the magnetic field at the center of the loop 
dotted into the distance (£), given by Eq (34) below, to the 
current carrying wire. 

The equation for the magnetic field of interest at a point 
(R,0,$) can be expressed as: 

B(R,0,®) = B( Ro , 6o ,<to ) + (del B)*£ (33) 


where 

£ = (R-Ro)R + (6-00)0 + ($-®o)3 (34) 


also 

B(R,0,«) s Bb(R,0,*)R + Bt(R,0,4»)0 + Bp(R,0,<&)4> (35) 

Expanding each scalar component: 

Br (R,0,*)=Br (Ro,0o,*0 )+3 Br/3R| oAR+3Br/S 0|o A6 +BBr/&$| oA$ (36) 
Bt (R,0,«)=Bt (Ro ,0o ,®o ) + dBT/aR|oAR+aBT/S0|oA0+aBT/0<I>|oA<I> (37) 
Bp (R, 6,4) =Bp (Ro ,0o ,4»o ) + 3Bp/3R} o AR+3 Bp/S 0 | o A6 +SBp| o A® ( 38) 
where AR = R-Ro and A0 = 0-0o . 

However Bp =0, which eliminates those four terms from 
consideration. Also Bb and Bt do not contain a phi term 
which eliminates two additional terms. Finally the terms for 
the field at the center of the loop will integrate to a zero 


a 
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net force—since they represent the loop in a uniform 
magnetic field—and need not be considered. 

The remaining radial and theta components of the magnetic 
field and their partial derivatives are expressed and 
evaluated as follows: 


&B*/0R = 0/3R[-2(Re/R) 3 g° cos( 0 ) ] 
3B*/00 = 3/ae[-2(Re/R) 3 g®cos(6)] 
BBt/BR = d/BR[ -(Re/R ) 3 g° sin( 0 ) ] 
8Bt/B6 = 3/00 [-(Re/R) 3 g° sin(0)] 
where 


6Re 3 g° cos ( 0 ) /R 4 (39) 

2(Re/R) 3 g° sin(0) (40) 

3Re 3 g° sin( 0 ) /R 4 (41) 

-(Re/R) 3 g° cos(0 ) (42) 


Re = earth’s radius 
R = orbit radius 
g° = first Gaussian coefficient 
0 = colatitude 

The evaluation of the gradient still requires the 
determination of AR and A0. From Figures 13-16 these can be 
seen to be ±rsin(Q)sin(defined orientation angle) and 
±(r/R)sin(Q)cos(defined orientation angle) respectively. In 
addition dL can be seen to be rdQ. 
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Figure 13. Case I: Current Loop and 

Orientation Relative to the 
Geomagnetic Field and Definitions 
of dL» AR and A6. 



Figure 14. Case II: Current Loop and 

Orientation Relative to the 
Geomagnetic Field and Definitions 
of dL• AR and A6. 
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Figure 15. Case III: Current Loop and Orientation Relative 
to the Geomagnetic Field and Definitions of dL, 
AR and A0. 



Figure 16. Case IV: Current Loop and Orientation Relative 
to the Geomagnetic Field and Definitions for dl*, 
ARi and A6. 
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When the current loop is aligned with its dipole parallel 
to the radial component of the geomagnetic field, the current 
loop’s orientation is approximately parallel to the earth’s 
surface. The periphery of the loop is assumed to be at the 
same distance as the center of the loop from the center of 
the earth. At a point exactly over the geomagnetic pole the 
radial component of the geomagnetic field yields a force that 
is directed radially with respect to the current loop and is 
therefore self-cancelling. Also, exactly over the magnetic 
pole the Bt component is zero. At equal, small distances 
away in all directions parallel to the geomagnetic equatorial 
plane the magnitude of Bt is small and equal at all points on 
the loop. The Bt component of the geomagnetic fields 
provides a force in the positive or negative radial direction 
(depending on the direction of current flow) with respect to 
the center of the earth. The Bt component at a location 
exactly centered over the pole therefore provides a net 
effective force of 2ariBT acting on the satellite. As the 
satellite moves away from the pole the net effective force 
will have contributions from both components of the 
geomagnetic field. 

The following derivation of the components of the force on 
the current loop makes the assumption that for each djt of the 
current loop the component of the geomagnetic field is 
parallel but of sightly different magnitude. Exactly over 
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the pole and for distances close to the pole this is not 
precisely correct. However, the convergence of components of 
B at the extreme points of the loop (1000 meter radius 
current loop in a 600 kilometer orbit) is less than 0.5 
degrees at one degree from the pole and less than 5 degrees 
when 0.1 degree from the pole. Figures 12 and 13 define the 
components of the geomagnetic field, orientation of the 
current loop and axes of the coordinate systems for the 
current loop and the geomagnetic field. 

The force on each infinitesimal of the current loop is 

d£ = idL x fi (43 ) 

The vector dl± can take one of two forms depending on the 
direction of current flow in the loop. For purposes of this 
derivation the current will be considered in a clockwise 
direction when the satellite is viewed from the earth’s 
center. 

The infinitesimal force expressed in determinent form is then 

d£ = i dL R 0 $ 

0 -cos(Q) -sin(Q) 

Br Bt 0 

Therefore, 

d£ = idL[B t sin(Q)R - BRsin(Q)8 + Brcos(2)$] (44) 

Each component of the force is given by the respective 
component in the preceding equation. The force in the radial 


direction (Fr) is given by: 

F* =JdFnR = Jirsin(Q)dQBr 


(45) 
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Substituting Eq (37) for Bt, and noting from Figure 13 that 


AR = 0 


( 46 ) 


and 


A6 = -(r/R)sin(2) (47) 

the integral becomes: 

J ir[-(Re/R) 3 g°cos(6)](-r/R)sin(2) sin(2)d2 (48) 

R and 9 define the location of the loop and i, r, g°, and Re 
are constants. Consolidating and pulling the constants out 
of the integral yields: 

ir 2 Re 3 g° cos ( 0 )/R 4 ^sin 2 ( 2 )dQ (49) 

When the integral is evaluated around the loop ( 2u) the force 
in the radial direction is: 


Far = iitr 2 Re 3 g°cos( 0 )/R 4 


(50) 


The force in the theta direction (Ft) is given by: 

Ft =Jd¥tQ = J-irsin( 2 )d2BR (51) 

Substituting Eq (36) for Br and recalling Eq (46) and (47) 
the integral becomes: 

y^-ir2(Re/R) 3 g°sin(9)(-r/R)sin(Q)sin(2)d2 (52) 

Pulling the constants out of the integral yields: 

2ir 2 Re 3 g°sin(9)/R 4 ^sin 2 (2)d2 (53) 

When evaluated around the loop the force in the theta 
direction is: 


Ft i = 2 ( inr 2 g° Re 3 /R 4 ) sin( 8 ) 


(54) 
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The force in phi direction (Fp) should be zero and the 
following evaluation confirms this. The force is represented 
by: 

Fp = fdFr* = J ircos ( Q )dQBn (55) 

Substituting Eg (36) for Br and recalling Eq (46) and (47) 
the integral becomes: 

Jir2( Re/R) 3 g° sin(6)(-r/R )cos( 2 )sin( 2 )dS2 (56) 

Consolidating terms and pulling constants to the left of the 
integral yields: 

-(2ir 2 g°Re 3 /R 4 )sin(6 )Jcos(2)sin(2)d2 (57) 

The integral, when evaluated around the loop (2it) yields zero, 
leaving a zero force in the phi direction. 

The total force is then: 

( 58 ) 

The direction of the force is dependent on the direction 
of the current flow. When the current is clockwise as seen 
from the center of the earth, then the radial, theta, and 
total forces vary with co-latitude as shown in Figure 17. 

The magnitudes are representative of a current loop carrying 
10 amps and a radius of 1000 meters. 
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Case II: Loop Perpendicu lar to Earth’s Surface 


With the current loop aligned with its dipole parallel to 
the theta component of the geomagnetic field (Bt), the 
current loop will maintain an orientation perpendicular to 
the earth’s surface as shown in Figure 12. Directly over the 
pole (Bt =0) the net force can be directed entirely in the 
direction of orbital motion by proper selection of the 
direction of current flow. As the satellite moves away from 
the pole, components of force are added in the radial 
direction. Again, for the following derivation, the 
assumption is made that the Bi or Bt components at the 
extremities of the current loop are parallel. Figures 12 and 
14 define for Case II the components of the geomagnetic 
field, orientation of the current loop, and axes of the 
coordinate systems for the current loop and the geomagnetic 
field. 

The force on the current loop is again given by Eq (43). 
However, since the loop now lies in the R-4> plane, dX. is given 
by dL[cos(Q)R 00 -sin(Q)$]. Expanding the cross product, the 

infinitesimal force is 

d£ = i dL [BTsin(Q)R - Bssin(Q)0 + Btcos(Q)$] (59) 

The procedure for determining the force in the radial and 
theta directions is the same as shown in Case I. 
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Substituting Eq (37) for Bt and noting from Figure 14 that 

A0 = 0 (60) 

and 

AR = rsin(Q) (61) 

the force in the radial direction is: 

(62) 

Which is zero over the pole (9=0) and a maximum at the 
equator as expected. 

Substituting Eq (36) for Br and recalling Eq (60) and (61) 
the force in the theta direction is: 

(63) 

The force in the phi direction is again zero because the 
integral ultimately contains the sin(Q)cos(Q) term. 

The total force acting on this orientation is then 

(64 ) 

Figure 18 shows the forces for Case II as they vary with co¬ 
latitude. The magnitudes are again representative of a 
current loop carrying 10 amps with a 1000 meter radius. 
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Figure 18. Magnetic Forces for Case II (600 Kilometer Orbit) 












For this third case the current loop is constantly 
changing its orientation with respect to the earth’s surface 
so as to maintain a constant orientation relative to the 
instantaneous geomagnetic field vector. The added complexity 
of this case is offset by the merit this orientation has in 
minimizing the torque experienced by the current loop. As 
can be seen in Figure 12 the satellite will make two complete 
rotations about its center of mass for each orbit about the 
earth. Identical to Case I when directly over the pole, it 
becomes Case II when at the equator. 

The determination of the force again assumes that the 
component of the geomagnetic field affecting each section of 
the loop is parallel for every section of the loop. The 
angle beta (0) defines the orientation of the loop at any 
given point in its orbit and, as can be seen in Figure 12 and 
15, is the angle between the current loop’s magnetic dipole 
(or instantaneous geomagnetic field vector) and the radius 
vector to the current loop from the center of the earth. 

The force on the loop is again given by Eq (43) and the 
expansion of the cross product is: 

dJE = i dL{BTsin(Q)R - BBsin(Q)9 

+ [ Bt cos ( Q ) sin( 0 ) + Br cos ( Q )cos ( 0 ) ]<&} (65) 

The force in the radial direction is given by Eq (45) 

Fb = J dFsR = ^irsin( Q )dQBT (45) 
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Substituting Eq (37) for Bi and noting from Figure 15 that 


AR = rsin(Q)sin(0) 


(66) 


and 


A6 = -(r/R)sin(Q)cos(0) 


( 67 ) 


the integral becomes: 

^irsin(Q)dfi{3(Re 3 /R 4 )g°sin(0)[rsin(Q)sin(0)] 

- (Re/R) 3 g°cos(9)(-r/R)sin(Q)cos(0)} (68) 

Forming two integrals, consolidating terms and pulling the 
constants out of the integrand yields: 

3sin(9)sin(0)(ir 2 g° Re 3 /R 4 )Jsin 2 ( Q )dQ 

+ cos(9)cos(0)( ir 2 g°Re 3 /R 4 )jsin 2 (C )dQ (69) 
Evaluating the integrals around the current loop yields a 
force in the radial direction of 


Fb 111 =(inr 2 g° Re 3 /R 4 )[3sin(0)sin(0)+cos(0)cos( 0)] 


( 70 ) 


Which gives the same results as Case I over the pole and Case 
II over the equator. 

The force in the theta direction is again given by Eq 
(51): 

Ft = f dF T 0 ~j- irsin(Q)dQBR (51) 

Substituting Eq (36) for Br and recalling Eq (66) and (67) 
the integral becomes: 

J- irsin(2)dQ{6(Re 3 /R 4 )g°cos(0)(rsin(0)sin(01] 

+2(Re/R) 3 g°sin(0)[(-r/R)sin(2)cos(0)]} (71) 





Forming two integrals, consolidating terms and pulling the 
constants outside the integrals yields: 

-6cos(0)sin(0)(ir 2 g° Re 3 /R 4 ) Jain 2 (Q)dQ 

+2sin{ 0)cos(0)(ir 2 g° Re 3 /R 4 ) fsin 2 (Q)dQ (72i 
Evaluating the integrals around the current loop gives the 
force in the theta direction as 


Ft 111 =2(inr 2 g° Re 3 /R 4 ) [sin(0)cos(0)-3cos(0)sin(0)] (73) 


Which again gives the same result as Case I over the pole and 
Case II over the equator. 

The force in the phi (<6) direction again equals zero 
because each integral contains a sin(Q)cos(2) term. The total 
force is given by 


( 74 ) 

The radial, theta, and total force for Case III are shown 
in Figure 19 for the same representative current loop as in 
Case I. 
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Figure 19. Magnetic Forces for Case III (600 Kilometer 
Orbit) 
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Case IV:_An laertially Oriented Current Lpqp 


In this situation the plane of the loop is maintained 
parallel to some inertial orientation. As long as the 
orientation can be defined in terms of the radial and theta 
components of the geomagnetic field lines, a solution can be 
found. For the purposes of illustration for this thesis, the 
orientation will be maintained parallel to the equatorial 
geomagnetic plane. Gamma (T ) is the angle between the plane 
of the loop and the radial direction and will be used to 
define the orientation as shown in Figure 12. When over the 
pole the orientation is the same as Case I. At the equator 
it resembles Case II or III but with the current flowing in 
the opposite direction. 

The force is again given by Eq (43). Gamma is 90 degrees 
over the pole and zero at the equator. Therefore T = 90-6 so 
cos(D = cos(9O-0) = sin(6) (75) 

sin(D = sin(9O-0) = cos(0) (76) 

The expansion of the cross product matrix can therefore be 
expressed 

d£ = idL{BTsin(Q)R-B h sin(Q)6 

+ [-BTCOs(Q)sin(0)+Bjicos(Q)cos(0)]<&} (77) 
The force in the radial direction is given by Eq (45) 

F* = fdF*R - J irsin(Q)d£2 Bt (45) 
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Substituting Eq (37) for Bt and noting from Figure 16 that 

AR = -rsin( 2 )cos( T) 

= -rsin(2)cos(90-6) 

= -rsin(2)sin(6) (78) 

and 

A0 = (-r/R)sin(2)sin(D 

= (-r/R)sin(2)sin(90-9) 

= (-r/R)sin(2)cos(9) (79) 

the integral becomes: 

Fr = irsin(2)d2{3(Re 3 /R 4 )g°sin(0)[-rsin(2)sin( 0)] 

+(-Re/R) 3 g°cos(0)[{-r/R)sin(2)cos(0)]} (80) 
Resolving into two integrals, combining terms and pulling 
constants out of the integrand yields: 

Fh = -3sin 2 ( 0 ) ( ir 2 g° Re 3 /R 4 )^sin 2 ( 2 )dQ 

+ cos 2 (0)( ir 2 g°Re 3 /R 4 )Jsin 2 (Q)d2 (81) 

Evaluating the integral around the current loop yields the 
force in the radial direction as 

(82 ) 

The force in the theta direction is given by Eq (51): 

Ft = JdFt 9 = J- irsin(Q)d2BR (51) 

Substituting Eq (36) for Br and recalling Eq (78) and (79) 
the integral expands to 

Ft = ^-irsin(2)d2{6(Re 3 /R 4 )g°cos(0)[-rsin(2)sin(9)] 

+ 2(Re/R) 3 g° sin(0)[(r/R)sin(2)cos(0)]} (83) 








T 


Resolving into two integrals, combining terms, and pulling 
constants out of the integrand yields: 

Ft = 6cos( 6 )sin(6 ) ir z g° Re 3 /R*Jsin 2 (2 )d2 

+ 2sin( 0 )cos ( 0 ) ir 2 g° Re 3 /R 4 ^sin 2 (Q )d2 (84) 

Evaluating the integral around the loop gives the force in 
the theta direction 


Ft iv = 8cos(0)sin(0)(inr 2 g°Re 3 /R 4 ) 


(85) 


The force in the phi (<t) direction is zero due to a sinlQ)cos(Q) 
term in the integrand. The radial, theta and total forces 
are shown in Figure 20 for Case IV. The magnitude is again 
representative of a 1000 meter radius current loop carrying 
10 amps. 
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V. APPLICATION 

The force alone tells little about the effect the current 
loop and geomagnetic field interaction have on a satellite. 
The force this technique can apply to a space system must be 
compared to other space propulsion systems. 

Various characteristics are used to compare thruster 
performance. The generally recognized measures are: specific 
impulse (Isp) - an overall performance comparison standard, 
characteristic exhaust velocity (c* ) - a measure of 
propellant combination and combustion chamber design, and 
thrust coefficient (Cr) - a measure of nozzle performance 
(Sutton, 1986:50,55). Since action-at-a-distance systems 
have no exhaust products and an infinite specific impulse, 
these traditional comparisons have little meaning. Some 
other comparison must be established. Other possibilities 
include total impulse (jFdt), power-plant specific mass 
(kilograms per kilowatt), or acceleration capability (Sutton, 
1986:21, Corliss, 1960:4). 

From the standpoint of acceleration capability, the 
obvious comparison is with other low-thrust systems such as 
ion propulsion, however a comparison with chemical rockets 
can be made based on total (lifetime) impulse. Naturally, 
comparisons will have to be made for operation in similar 
environments - for instance "gravity-free", drag-free space. 

Practically speaking it is the acceleration this force can 
produce to perturb the satellite that is of interest. The 
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acceleration stems from the relationship given by Newtons 
second law, that is, F = ma . Therefore, satellite system 
mass (m) must be known to assess the acceleration (a) 
available from a given force (F). 

Some baseline satellite characteristics will be 
established. Assuming a satellite operating payload of 
approximately 3000 kilograms, the mass of the power supply, 
the current loop, the support structure, and the power 
conditioner must be added. The current loop mass is a 
function of the density of its material and size. Comparing 
nine gauge (2.91 millimeter diameter) wire for two common 
conductors, copper and aluminum, the following 

characteristics are given {CRC, 1973:F137-143): Copper has a 
current capacity of 30 amperes, a resistivity of 2.60 ohms 
per kilometer (at 20'C), and a mass of 58.98 kilograms per 
kilometer. Aluminum has a current capacity of 25 amperes, a 
resistivity of 4.26 ohms per kilometer, and a mass of 17.90 
kilograms per kilometer. A 1000 meter radius loop of copper 
would have a mass of 371 kilograms and the same loop in 
aluminum would be 113 kilograms. From the standpoint of cost 
to launch, lighter materials are preferable to heavier. The 
high currents required in this system require thicker wire to 
reduce resistance. The more conductive wire will reduce 
ohmic heating, improving efficiency. A measure of merit used 
by Prall (1987:38) would be to look for the largest ratio of 
conductivity to density. It is apparent that little is lost 
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in terms of resistance and current capacity if aluminum wire 
is used rather than the heavier copper. 

The power generation could be from a nuclear reactor or 
solar panels. Assuming solar cells were used, a reasonable 
attainable goal for solar panels is 250 watts per square 
meter and 20 kilograms per kilowatt (Sutton, 1986:498). To 
generate 25 kilowatts of power for the satellite systems and 
to drive a current of ten amps in the loop would require 100 
square meters of solar cells and 500 kilograms of additional 
mass. Since optimal alignment with respect to the sun may 
not always be maintained and degradation of the solar cells 
will occur over time, doubling this to 200 square meters and 
1000 kilograms provides a margin for continuing capability. 

A support structure for the current loop, a power 
conditioner, and accessories for the satellite could mass 
5885 kilograms for a system mass of 10,000 kilograms. Using 
representative values of ten amperes for the current (i) and 
1000 meters for the current loop radius (r) and the satellite 
mass arrived at above, the acceleration capability ranges 
from zero to 6x10" 4 g’s and a watts per pound ratio of 
approximately one. Point A on Figure 21 compares this system 
to other typical values for different types of rocket 
engines. 

Another comparison can also be made using total impulse. 
Total impulse (It) is the thrust force (F) integrated over 
the operating time (t); It = /Vdt , if the thrust is 
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constant and there are no stop and start transients, It = Ft. 
From Figure 17-20 an average magnetic interactive force of 
two (2) newtons appears reasonable and assuming a twenty year 
lifetime, the total impulse is 1.26x10®Kg-m/sec. 



Figure 21. Operating Regimes of Various Thruster Systems 
(Adapted from Sutton:32) 


Using approximate values for thrust and burn time for 
various chemical rockets from Sutton 

(1986:31,196,210,216,263) and electric thrusters from Fearn 
(1982:160) total impulse can be calculated and is shown in 
Table I. 
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Table I. Comparison of Thrust, Burn Time, and Total Impulse 
for Various Rocket Engines 

Rocket Engine Thrust (Nt) Sum Tipe ( 9 ) To.tal Impulse 


Saturn F-l 

7.34E6 

165 

1.21E9 

Saturn J-2 

1.02E6 

500 

5.12E8 

Shuttle Main 

1.89E6 

480 

9.07E8 

Minuteman 

8.90E5 

55 

4.89E7 

Lunar Module 

4.63E4 

800 

3.70E7 

General RCS 

6.67E1 

3600 

2.40E5 

Shuttle RCS 

1.11E2 

8.2E4 

9.10E6 

Kaufman 8 cm 

4.90E-3 

5.4E7 

2.65E5 

Kaufman 30 cm 

2.10E-1 

3.6E7 

7.56E6 


By this standard of comparison the magnetic field 
interaction appears very favorable, however, some of the 
comparisons are unrealistic in that the first five listed are 
primary propulsion designed to boost a vehicle out of a 
gravity well—something magnetic field interaction simply 
cannot accomplish. A more appropriate comparison would be to 
look at impulse deliverable over one orbit. Given an 
initially circular orbit of 600 kilometers, the orbital 
period is 5800 seconds. Table II shows impulse delivered in 
5800 seconds for several thrusters. 


Table II. Impulse Over One Orbit for Various Thrusters 


Thruster 
Magnetic 
Generic RCS 
Shuttle RCS 
Kaufman 8 cm 
Kaufman 30 cm 


Thrust (Nt) 
2.0 
66.7 
111.2 
0.0049 
0.210 


Total Impulse 
1.16E4 
3.87E5 
6.44E5 
2.84E1 
1.22E3 
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An external acceleration, that is, an acceleration other 
than the acceleration due to the central gravitational field 
that defines the orbit, will cause the parameters that define 
the orbit to vary with time. Ideally all disturbing forces 
must be determined and the net acceleration vectors 
calculated to specify how the orbit will change in the 
presence of those forces. Some of these disturbing forces 
that act on an orbiting satellite include aerodynamic drag, 
solar radiation pressure, earth’s oblateness, induced 
electromagnetic force, solar wind, cosmic dust etc. 

For this analysis, a comparison will be made between the 
magnetic force acting on the satellite and aerodynamic drag 
and solar radiation pressure. An estimate can be made of the 
change in orbital radius by equating the work done by the 
force over some portion of the orbit to the change in the 
specific mechanical energy of the orbit. In addition, 
Appendix A shows how the magnetic force can be applied using 
the Variation of Parameters technique to determine how the 
orbit can be altered by the magnetic force. 

Magnetic Force vs Aerodynamic Drag 

Aerodynamic drag forces are tangent to the orbit and 
opposite to the velocity. Since atmospheric particle density 
is a function of altitude, the aerodynamic drag forces are 
most noticeable at satellite perigee. Qualitatively, drag 





effects are similar to an in-plane transfer maneuver applied 
at perigee: apogee height will be reduced as well as the 
semi-major axis and eccentricity. Perigee height, argument 
of perigee and inclination will remain approximately the 
same. The latter does not change since the aerodynamic drag 
effects are applied in-plane. 

If atmospheric rotation is ignored—which is consistent 
with assuming a non rotating earth--the speed of the 
satellite relative to the atmosphere equals the inertial 
velocity. The aerodynamic drag force (Fa) can then be 
expressed as 

Fd = iCapv 2 A ( 86 ) 

where 

Cd = Drag coefficient 
p - Atmospheric density 
v = velocity of satellite 
A = frontal area of satellite 

There are various frontal area configurations that could be 

considered to determine the drag coefficient. The solar 

panels could face into the direction of motion; the current 

loop support structure could be edge on or perpendicular to 

the orbit motion etc. The drag coefficient is approximately 

one ( 1 ) when the mean free path of the atmospheric molecules 

is small compared to the size of the satellite and two ( 2 ) 

when the mean free path is large compared to the size of the 

satellite (Roy, 1965:230). NASA SP 8058 gives a conservative 

drag coefficient of 2.6 for flat plate shapes. Assuming the 

drag coefficient to be a very conservative 3.0, the frontal 
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cross sectional area of the satellite system to be a minimum 
of 80 square meters and a maximum of 300 square meters, Table 
III shows the estimated aerodynamic force from 200 to 1000 
kilometers. 

Table III. Aerodynamic Drag Force at Various Altitudes 
Altitude (Km) Density( Kg/m 3 ) Velocitvt m/s ) Force (Nt 18,0/300m 2 


200 

2.541E-10 

7.789E3 

1.603/6.0123 

400 

2.803E-12 

7.673E3 

1.716E-2/6.436E-2 

600 

1.137E-13 

7.562E3 

6.762E-4/2.536E-3 

800 

1.136E-14 

7.456E3 

6.568E-5/2.463E-4 

1000 

3.561E-15 

7.354E3 

2.003E-5/7.5UE-5 


Figure 22 shows how the total magnetic force varies with 
altitude for Case III. Again magnitudes are representative 
of a 1000 meter radius current loop carrying 10 amps. As can 
be seen, the magnitudes of the magnetic force are all greater 
than the aerodynamic drag at altitudes over 400 kilometers. 
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Solar Bad iali aa Ersaaitce 


Sutton (1986:117) states that the solar radiation pressure 
on a given surface of a satellite near the earth’s solar 
distance is 

p s 4.5E-6 cos(0)I(l-k«)cos(0) + 0.67kd] (87) 

where 

ks= specular coefficient of reflectivity 
kd= diffuse coefficient of reflectivity 

0 = angle between incident vector and surface normal vector 
Typical values of 0.9 and 0.5 for ks and kd respectively 
for the body of the satellite and 0.25 and 0.01 for the solar 
panels are given in Sutton and NASA SP 8027. Using an 
approximate satellite body area of 25 square meters and 
effective solar array area of 80 to 300 square meters the 
solar radiation pressure is 3.214E-4 to 1.0704E-3 newtons. 

This force begins to dominates aerodynamic drag forces 
between 600 and 800 kilometers as shown in Figure 23. 
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Applying low, continuously applied tangential thrust in 
the direction of satellite motion, potential energy rather 
than kinetic energy is added to the satellite. The satellite 
will progress to a higher, more energetic orbit, although one 
with less tangential velocity (Corliss, 1960:38). The work 
done by the force is the product of the force and the 
distance over which it is applied [FRd9] . This can be added 
to the specific mechanical energy of the orbit or equated to 
the change in energy of the orbit [d/dR(-p/2R) = p/2R 2 ]. 

Given that the satellite starts out in a circular orbit and a 
very low tangential force is applied, the satellite will 
follow a spiral trajectory to a higher orbit. The orbit 
will, however, continue to approximate a circle. Corliss 
gives the increase in orbit radius by the application of a 
small tangentially applied force as 

dR = 2FR« /go m( R/R* ) 3 d6 (88) 

where 

R = orbital radius (m) 

R*= radius of the earth (m) 

F = force (Nt) 

m = mass of the satellite (kg) 
go= acceleration of gravity (m/s 2 ) 

0 = angle along the orbit (radians) 

Introducing time from Vs = Rd0/dt, we have 

dt = (Rs/go )WZ (R/R. ) 3 ' 2 d0 (89) 
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Combining the two equations then gives 

t = gomR./F(go )W2 [l/(R 0 )W2 _ l/(R)l/2) (90) 

where 

R = initial radius of satellite’s orbit (m) 

Ro= final radius of satellite’s orbit (m) 

Corliss shows that a satellite with a thrust to weight 
ratio of 10" 4 (initially in a 500 kilometer circular orbit) 
can reach escape velocity in 92 days and 400 orbits. This is 
the same order of magnitude as the acceleration available due 
to the magnetic field interactive force. However, since the 
magnetic field strength drops off as the cube of the radius 
of the orbit, the applied force also decreases. While escape 
velocity cannot be achieved in ninety days, the orbit radius 
can be changed by approximately 1000 meters over one half of 
an orbit for a ten amp, 1000 meter radius current loop, 
starting in a 600 kilometer circular orbit. This orbital 
radius change assumes a constant force of two newtons applied 
in the tangential direction. 

For the four orientations discussed in Chapter IV, the 
force varies in magnitude and is not totally applied in a 
tangential direction, therefore the orbit will not actually 
be a spiral but osculatory as it gains energy. If the 
orientation and current direction can be changed at will, 
higher average force levels are available. A more precise 
determination of the change in the orbital elements requires 
an extensive analysis using the variation of parameters. The 
appropriate equations are given in Appendix A. 






VI* .SfiNCLU S IPN 


Summary 

Several models exist to determine the geomagnetic field at 
any point around the earth. The simplest of these, and the 
one used herein to demonstrate the thesis concept, is the 
axial dipole. It was shown that the force acting on the 
current loop by the geomagnetic field can be calculated for 
any orientation of the current loop in orbit about the earth. 

Given the assumptions in Chapter I and the constraints 
they represent, a number of conclusions can be drawn. The 
first, that a small net force exists on a current loop that 
extends over some gradient of the geomagnetic field, was 
anticipated, as well as the fact that it would scale 
proportionally to the area of the loop and the current the 
loop would carry. The fact that the net force would vary to 
the extent that it does with position and orientation was not 
expected. 

As was shown in Table III and Figure 22, the magnetic 
force exceeds aerodynamic drag above approximately 300 
kilometers. Even as solar radiation pressure begins to 
exceed aerodynamic drag above 700 kilometers, it, in turn, is 
still less than the magnetic force available for most of an 
orbit up to at least 1600 kilometers. The preeminent 
conclusion to this thesis, therefore, is that a current loop 
can be applied as a space thruster to change the specific 
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energy of a satellite and hence the satellite’s orbital 
elements. 

While several simplifying assumptions were made, they did 
not detract from the proof-of-concept. For instance, the 
simplification of the earth’s magnetic field to a simple 
axial dipole was made for ease of analysis, but the method of 
analysis would be the same if an octopole or greater 
expansion of the infinite series was used. With only the 
first term of the geomagnetic field expansion being used 
there was no phi component of the magnetic field. 

Furthermore, only polar orbits were considered since there 
would be no force component in the phi direction, however, 
any orbit can be considered. The four orientations 
represented a building block approach to the analysis, but 
their selection covers all possible orientations and provide 
a stepping off point for the inclusion of the effects of 
torque on the loop. The non-rotating earth simplification was 
made for the sake of atmospheric drag calculations. 

Some of the assumptions were made due to a lack of 
definitive engineering data and experience in the academic 
and engineering worlds. Specifically, the assumption of 
rigidity and power requirements and the dynamics of these 
structures. For proof of concept, however, these assumptions 
are tolerable, though subject to possible debate. The 
dynamics of large scale space structures and power generation 
in space will be answered in greater detail as the nation 
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ventures into the space station era of space exploitation and 
exploration. 


Recommendations 

Should this type of action-at-a-distance thruster actually 
be employed, a computer generated solution should carry the 
geomagnetic model to an octopole or even include all 80 terms 
of the IGRF model with coefficients updated by a follow-on to 
MAGSAT. This would be necessary for an accurate 
determination of the force the current loop would experience, 
as well as account for the changes that the geomagnetic field 
undergoes as opposed to a simple axial magnetic dipole. 

Since one assumption made was the ability to maintain any 
desired orientation throughout the orbit, additional work 
should be done to determine the net force as the loop rotated 
due to the torque it would experience in the geomagnetic 
field. This would probably be oscillatory and approach a 
stabilization orientation similar to Case III. Another 
useful analysis would be to solve for the orientation that 
maintains the greatest force for any given position. 

Since the effective use of this type of action-at-a- 
distance force depends on the advent of large scale space 
structures and large amounts of power, an initial application 
could perhaps be considered for an orbiting power generation 


satellite. 





Variation of Parameters is a computational method for 
calculating orbital elements. In the presence of a small 
perturbative force the orbital elements will change slowly with 
time. While other methods of calculating changes to an orbit 
maintain a constant reference orbit until rectification occurs, 
with variation of parameters the orbit changes are calculated 
continuously. 

The following steps and equations (adapted from Bate and 
others, 1971:396-407,72) outline the procedures necessary to 
determine the actual change that will take place in a 
satellite’s orbit when subject to a perturbative force. 

1. Set initial conditions at time t = to for position (jr) 
and velocity (v) in the (I,J,K) reference frame (Figure 

11 ). 

2. Calculate the fundamental vectors (h,n>e) at time t = to 
where 

& = the angular momentum (h = £ x v) 

B = the node vector (n = K x h) 

§. = the eccentricity vector (e = 1/p[(v 2 -p/r)£ - 

( r • v )v ] ) 

3. Calculate the orbital elements at time t = to 
where 

a = the semi-major axis 


e = the eccentricity 






i = the inclination 


Q = longitude of ascending node 
w = argument of periapsis 
These are shown in Figure 11. 

4. Compute the perturbative force at time t = to and 
transform to the (R,S,W) system. These force components 
are given in Eq (50), (54), (62), (63), (70), (73), 

(82), and (85). 

5. Compute the right hand side of da/dt, de/dt, di/dt, 
dQ/dt, dw/dt, dMo/dt (as shown in steps 5-1 through 5-7 
below). 

6. Numerically integrate each equation from step five one 
time increment and step time counter. 

7. Add the changes from step six to the old elements from 
step three. 

8. Calculate the new position and velocity vectors 
(jE and v) from the new elements. 

9. Go to step four until reaching final time step. 

Each of the preceding nine steps is expanded in the following 
series of steps and substeps: 


Set initial 

conditions 

at t = to 


£ = ril + 

A A 

r j J + rs K 

r = ( ri 2 

+ rj 2 + rn 2 ) 1/2 

A 

A A 


+ VJ 2 + VK 2 ) 1 ^ 2 

V = VI I + 

VJ J + VJ J 

v = (vi 2 
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2. Calculate the three fundamental vectors 

2-1. h = r x v = I J K 

ri rj ric 

VI vj vn 

A 

= (rj vk - vjrn)I - (ri vu - vitk)J + ( ri vj - vi rj ) 

= hi I + hj J + hKK h = (hi 2 +hj 2 + hit 2 ) 1 / 2 

2-2. n = k x h = I J K 

0 0 1 

hi hj hs 

= -hjl + hi J 

= m I + nj J + nit K n = (m 2 +nj 2 + me 2 ) 1 ^ 2 

2- 3. e = l/p{[v 2 - (p/r)]r - (r • v)v]} 

= l/p{[v 2 - (u/r)](ril + r j J + ricK) 

-[ (rivi + n vj + tkvk)(vi + vj + vk ) 1 

where p = GM 

= ei I + ej J + eic K e=(ei 2 +ej 2 +eit 2 ) 1/2 

3. Calculate the orbital elements 

3- 1. p(semi-latus rectum) = h 2 /p 

3-2. e = |a| 

3-3. a = (1-e 2 )/p 

3-4. i = arccos hK/h 

3-5. Q = arccos m/n (if nj>0 ,then 2 <180') 

3-6. w = arccos (a • s)/ne =(m ei I + nj ej J + n« eiKl/ne 

= (-hj ei I + hiejJ)/ne (if es >0 , then w<180°) 

3-7. u 9 (true anomaly at epoch) = (s. * l)/er 

A A A 

= (eiril +ej rj J + ei nK)/er 

(if £*v>0 ,then <180°) 
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3-8. uo(argument of latitude at epoch) = (q • r)/nr 

A A a 

= (niril + nj rj J + nirKK)/nr 

= (-hjn + hirj)/nr (if rs >0 ,then uo<180°) 

3-9. lo(true longitude at epoch) = Q + w+k#=2 +uo 


4. Calculate the perturbative force 

4-1. K = (inr 2 g° Re 3 /R 4 ) 

4-2. FR1 = K cos(6) = F* 

4-3. FT1 = K 2sin(0 ) = Fs 

4-4. FR2 = K 3sin(0 ) = Fr 

4-5. FT2 =-K 6cos(0 ) = Fs 

4-6. FR3 = Kl3sin(0)sin(0) + cos(0)cos(0 ) ] = Fr 

4-7. FT3 = K[2sin(0)cos(0) - 6cos(0)sin(0)] = Fs 

4-8. FR4 = K( cos 2 (0) - 3sin 2 (0)] = Fr 

4- 9. FT4 = K[8cos(6)sin(8)] = Fs 

5. Compute the right hand side of the changes in the 
orbital elements. 

5- 1. da/dt = [2esin (v )F*/n(1-e 2 ) 112 ] + [2a(1-e 2 ) 11 2 Fs /nr ] 

5-2. de/dt = [(1-e 2 ) 1 > 2 sin (v )Fe/ na] 

+ [ (1-e 2 ) 1 f 2 /na 2 e ] { [a 2 ( 1-e 2 )/r ] - r}Fs 
5-3. di/dt = rFwcos(u)/na 2 (1-e 2 ) l < 2 ( = 0 since Fw= 0) 

5-4. dQ/dt = rFw sin(u )/na 2 (1-e 2 ) x > 2 sin', i ) (=0 since Fw=0) 

5-5a. (dw/dt)a = -(1-e 2 ) 1 / 2 cos( v) Fr/ nae 

5-5b. (dw/dt)s = p/eh{sin (v ) [1 + l/[l + ecos (v )]}Fs 

5-5c. (dw/dt)w = -rcot(i)sin(u)Fw/na 2 (1-e 2 ) 11 2 

5-6. dn/dt = (-3p/2na 4 )da/dt 
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5-7. dMo/dt= -l/na{<2r/a) - [(1-e 2 )/e]cos (v ) }Fr 

- [(1-e 2 )/nae]{1 + r/[a(1-e 2 )]}sin( v )Fs 

- t(dn/dt) 


6 . 


7. 


8 . 


8-1 

8-2 

8-3 

8-4 


R 


9. 


Numerically integrate each equation in step five above 
one time step. 

Add the changes arrived at in step six to each of the 
old elements from step three. 

Calculate new position and velocity vectors from the new 
elements arrived at in step seven. This is most easily 
done by calculating r and v in the perifocal coordinate 
and transforming back to the geocentric-equatorial 
(I,J,K) system. 

r = p/[1+ ecos( V )] 

A A 

j: * rcos(iOP + rsin(v)Q 
v = ( p/r ) 1 1 2 {-sin( v )P + [e + cos(l>)]Q] 

Transform to (I,J,K) 


ai 


ap 

aj 

=R 

a«j 

ax 


aw 


Let c = cos 
s = sin 

: c(Q)c(w)-s(Q)s(w)c(i ) -c(Q)s(w)-s(Q)c(w)c(i) s(fi)s(i) 
c (w) s( 2 )+c ( Q ) s( w )c( i ) -s(2)s(w)+c(2)c(w )c(i ) -c(2)s(i) 
s(w)s(i) c(w)s(i) c(i) 

Go to step four until reaching final orbit or last time 

increment. 


The steps are summarized in the following flow chart. 
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[ Figure 

i 


24. 


Variation of Paraaeters Flow Chart 
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Appendix B: Program for Calculating and 
Plotting Magnetic ..Forsg 


The force calculations were made through the use of a 
software package called Powerpack. Created as a statistical 
applications software it is also capable of performing various 
formula calculations and has a simple but useful plotting 
routine. The following is the core of the magnetic force 
computation program. 


•THIS PROGRAM CALCULATES THE GEOMAGNETIC FIELD AS A FUNCTION OF THE RADIAL; 
•DISTANCE FROM THE CENTER OF THE EARTH AND COLATITUDE USING A SIMPLE; 

•AXIAL DIPOLE MODEL FOR THE MAGNETIC FIELD. THIS PROGRAM THEN CALCULATES; 
•THE FORCE ON A CURRENT LOOP IN THE GEOMAGNETIC FIELD. THE CURRENT i IS; 
•ENTERED IN AMPS, THE CURRENT LOOP RADIUS r IS ENTERED IN METERS AND THE; 
•ALTITUDE h IS ENTERED IN METERS. THERE ARE FOUR ORIENTATIONS TO BE; 
•CONSIDERED AND TWO POSSIBLE CURRENT DIRECTIONS; 


•THIS FIRST SECTION DETERMINES WHAT VALUES ARE TO BE RETAINED UP TO A; 
•MAX OF 10, AND THE INITIAL INPUT PARAMETERS; 


RETAIN THETA(DEGREES),BETAd,KM,F3; 
prompt l,r,h,t; 
define Re * 0370000; 
define Rh * Re * h; 

--; 

•THIS SECTION CALCULATES THE GEOMAGNETIC FIELD AND ITS ORIENTATION; 


define Bt « -(0370000/Rh)*3«,3I*eln(t); 

define BR » -2*(0370000/Rh)*3«.31*ooe(t); 

define B - (Bt'2 + BR’2)'.S; 

define ■ * (Bt); 

define c = (BR); 

define Beta ■ ArcTan(e.c); 

•-; 

•THIS SECTION CALCULATES THE FORCE FOR THE DIFFERENT ORIENTATIONS; 


define K • (3 . 1 4 15fi»r*2 » 1 ». 3 1 »Re‘3/Rh‘4) ; 


•CASE 1: CURRENT LOOP PARALLEL TO EARTHS SURFACE; 


define FR1 = K»co«(t); 
define Ftl * K»2««in(t); 
define FI ■ (FR1"2 ♦ Ftl”2)‘.S; 


■CASE 2: CURRENT LOOP PERPENDICULAR TO EARTHS SURFACE; 

-- 

define FR2 ■ K*3»ain(t); 
define Ft2 « -K*fl*ooe(t); 

define F2 - (FR2*2 ♦ Ft2“2)*.S; ' 
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•CASS 3: CUBBEHT LOOP DIPOLE TAHOEHTIAL TO OEOMAOHETIC FIELD; 

* - 

define FB3 >Aba(K»(3*aln(t)«eln(Bet*)*coi(t)*eoe(Beta)); 
define Ft3 »Aba(K*(2*co#(Beta)*»ln(t)-8*cosCt)««ln(Beta)) 
define F3 « (FB3“2 ♦ Ft3*2)*.S; 

• - 

• CASE 4: CURREHT LOOP MAIHTAIXISO CONSTANT INERTIAL ATTITUDE 

• - - 

define FB4 ■ K*(eoe(t)»eoe(t)-3»aln(t)»aln(t)) 
define Ft4 * B»K*eoa(t)»eln(t) ; 
define F4 * (FB4*2 «• Ft4*2)*.3; 


•REVERSING THE CUBBEHT DIBECTIOH 


define FB1B > FB1»(-1); 

define FtlB * Ftl»(-1); 

define FB2B « FB2«(-1); 

define Ft2B « Ft2*(-1); 

define FB3B * FB3»(-1); 

define Ft3B * (-1)*ABS(2«K*(3*coa(t)»*in(Beta)-aln(t)*co*(Beta))) 

define FB4B « K»(3*cln(t)•eln(t)-coe(t)•eoe(t)); 

define Ft4B - Ft4*(-1); 

define THETA(DEOBEES) » (180«t)/3.1413827; 
define BETAd « 180 ♦ (180»Beta)/3.1415827; 
define KM * h/1000; 

LET HSMIN * 0; LET HSMAX * 380; 

LET VSMIN * 0; LET VSMAX « 4.3; 

LET LHGXMIN ■ 6; LET LEOYMIH >0.1; 

SYMBOL 0 7 30 42 31; 

YLABEL » FORCE (NEWTONS); 
run; 
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